Do Sweeping Effects Suppress Particle Dispersion in Synthetic Turbulence? 
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Synthetic models of Eulerian turbulence like so called "Kinematic Simulations" (KS) are often 
used as computational shortcuts for studying Lagrangian properties of turbulence. These models 
have been criticized by Thomson & Devenish (2005), who argued on physical grounds that sweep- 
ing decorrelation effects suppress pair dispersion in such models. We derive analytical results for 
Eulerian turbulence modeled by Gaussian random fields, in particular for the case with zero mean 
velocity. Our starting point is an exact integrodifferential equation for the particle pair separation 
distribution obtained from the Gaussian integration-by-parts identity. When memory times of par- 
ticle locations are short, a Markovian approximation leads to a Richardson-type diffusion model. 
We obtain a time-dependent pair diffusivity tensor of the form Kij(r,t) = Sij(r)r(r,t) where 5V, (r) 
is the structure- function tensor and r(r,t) is an effective correlation time of velocity increments. 
Crucially, this is found to be the minimum value of three times: the intrinsic turnover time T e dd y (r) 
at separation r, the overall evolution time t, and the sweeping time r/vo with vo the rms velocity. 
We study the diffusion model numerically by a Monte Carlo method. With moderate inertial-ranges 
like those achieved in current KS, our model is found to reproduce the t 9 ' 2 power-law for pair disper- 
sion predicted by Thomson & Devenish and observed in the KS. However, for much longer ranges, 
our model exhibits three distinct pair-dispersion laws in the inertial-range: a Batchelor ^-regime, 
followed by a Kraichnan-model-like t 1 diffusive regime, and then a r 6 regime. Finally, outside the 
inertial-range, there is another t 1 regime with particles undergoing independent Taylor diffusion. 
These scalings are exactly the same as those predicted by Thomson & Devenish for KS with large 
mean velocities, which we argue hold also for KS with zero mean velocity. Our results support the 
basic conclusion of Thomson & Devenish (2005) that sweeping effects make Lagrangian properties 
of KS fundamentally different from hydrodynamic turbulence for very extended inertial-ranges. 

PACS numbers: 47.27.Ak, 47.27.tb, 47.27.eb, 47.27.E- 



I. INTRODUCTION 

How particle pairs separate in a turbulent flow has 
been a central subject of turbulence research since the 
classical work of Richardson [l[ . Unfortunately, the phe- 
nomenon has proved quite difficult to investigate by nu- 
merical solution of the fluid equations and by controlled 
laboratory experiments, especially because of the very 
large Reynolds numbers required. Many studies have 
therefore employed "synthetic turbulence" or ensembles 
of random velocity fields with some of the scaling prop- 
erties of real turbulent velocities but which can be effi- 
ciently sampled even for very long scaling ranges. For 
example, papers have followed this approach and 

have reported substantial agreement of their numerical 
simulations with the predictions of Richardson, includ- 
ing the famous -law" for the growth in time of mean 
square pair separation distances. 

The validity of these results has been called into ques- 
tion, however. A paper of Chaves et al. Q pointed 
out that the use of synthetic turbulence to model Eu- 
lerian velocity statistics implies sweeping effects of large- 
scale eddies on particle motions that diverge with the 



Reynolds number. Those authors suggested to employ 
synthetic ensembles such as Gaussian random fields to 
model instead the turbulent statistics of Lagrangian ve- 
locities. In a simple one-dimensional Gaussian model 
of Eulerian velocities they found analytically that large- 
scale sweeping effects "localized" particle pairs and pre- 
vented them from separating. Subsequently, in a very 
interesting paper Q, Thomson & Devenish have pro- 
posed an intuitive picture how sweeping affects parti- 
cle dispersion in synthetic models of Eulerian turbu- 
lence. The key point is that large-scale eddies in real 
turbulence advect both particles and smaller scale ed- 
dies, while large-scale eddies in synthetic turbulence ad- 
vect only particle pairs and not smaller eddies. This fact 
implies that particle pairs at separations r in synthetic 
turbulence should experience rapidly changing relative 
velocities, as they are swept into new, statistically inde- 
pendent eddies. This occurs on a "sweeping" time-scale 
T S w( r ) ~ r /vo, where vq is the rms velocity set by the 
largest eddies in the synthetic ensemble. Thomson & 
Devenish assume a diffusion process of pair separations 
with an eddy-diffusivity K(r) ~ Su 2 (r)T sw (r) and Su 2 (r) 
the mean-square relative velocity at separation r. In an 
ensemble with Kolmogorov scaling Su 2 (r) ~ (er) 2 / 3 , this 
yields dr 2 jdt ~ K(r) ~ e 2 ^ 3 r 5 ^ 3 /vq and the solution 
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Note that this implies considerably slower growth than 
Richardson's f 3 -law [3(j. Thomson & Devenish argued 
for the above prediction in the case of a large mean sweep- 
ing, with vq replaced by the mean speed u. In the case 
of a zero-mean velocity ensemble, they argued instead 
for a t 9 / 2 -growth law, intermediate between t 3 and t 6 
(see our section 11111 below) . These predictions were sup- 
ported in Q by the numerical technique of "Kinematic 
Simulations" (KS) d-Q. The previous contrary numer- 
ical results were explained on various grounds, e.g. the 
use of an adaptive time-stepping scheme in Q which did 
not resolve the small sweeping time T sw (r) and its effect 
on particle dispersion. 

The issues raised by the paper of Thomson & De- 
venish have still not been fully resolved. The numerical 
simulations in Q used another form of adaptive time- 
stepping, which was suggested in Q to be responsible 
for the observation of a t 9 / 2 growth. Thomson & De- 
venish then repeated their simulations with a fixed small 
time-step and reported the same t 9 / 2 law [lj| ■ The most 



recent simulations of Nicolleau & Nowakowski for 
their longest scaling ranges show some evidence of the 
Thomson-Devenish sweeping effects, but the reported 
scaling laws are intermediate between those of Richard- 
son and of Thomson-Devenish and agree with neither 
theory. Thus, there is still considerable uncertainty in 
the literature regarding the validity of the Thomson- 
Devenish theory. The question is important, because 
synthetic turbulence is a useful testing ground for nu- 
merical and theoretical methods, and because compari- 
son of particle dispersion in synthetic and real turbulence 
illuminates the physical mechanisms of the latter. 

Because of the disagreement of the numerical simula- 
tions of different groups, it is useful to have analytic re- 
sults. The Thomson-Devenish arguments apply to a wide 
array of synthetic turbulence models, but Gaussian ve- 
locity ensembles are the most mathematically tractable. 
We therefore consider here the use of Gaussian random 
fields as models of turbulent Eulcrian velocities. More 
precisely, we take the advecting velocity field u(x, t) to 
be a Gaussian random field with mean u(x, t) and co- 
variance Cy(x, t; y, ,s) = (?4(x, t)u' 3 (y, s)) for the fluctu- 
ations u' = u — u. Specific models of interest are sim- 
ilar to those studied in 0, with u(x,i) = u indepen- 
dent of space and time and with covariance defined for 
0<a<2,0</3<2by 
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Here k\ = k 2 + l/L 2 and P„(k) = <% - kikj/k 2 is the 
projection onto the subspace of M. d orthogonal to k. The 
Gaussian random field u(x, t) so defined is statistically 
homogeneous in space, stationary in time, and solenoidal. 
The length L is proportional to the integral length-scale. 
The scaling properties of the model at scales smaller than 
L are similar to those of real turbulence. For example, 



the single-time covariance for r <C L is calculated to be 

C y (x,i;y,i) ~D L a 

-D^Hd + a-^-ahfjj+O^/L 2 ) (3) 

with r = x — y. See [12], p. 686. Kolmogorov 1941 dimen- 
sional scaling corresponds to the exponents a = /? = 2/3. 
We shall consider also in this paper Gaussian velocity 
models whose energy spectra coincide with KS models 
which have been studied numerically @-[ll| . The incom- 
prcssibility of these models will be used in an essential 
way, although much of our analysis applies to more gen- 
eral models, e.g. with any degree of compressibility. 

The principal results of this paper are as follows. For a 
general Gaussian model of Eulerian turbulence we care- 
fully derive the diffusion approximation for pair disper- 
sion assumed in the argument of Thomson-Devenish [8| , 
under the assumption of short memory times for parti- 
cle locations. We furthermore obtain a closed formula, 
cq. ([54"l) . for the 2-particle eddy-diffusivity in a general 
Gaussian model. For the specific models with covari- 
ance ^ we obtained more explicit results, which, under 
the conditions a < 1 and either /3 < 1 or frozen turbu- 
lence with D3 = 0, verify the Thomson-Devenish argu- 
ment about sweeping decorrelation effects. In particu- 
lar, we obtain under these conditions a 2-particle eddy- 
diffusivity tensor of the form Kij(r,t) = Sij(r)T(r,t), 
where (r) is the structure-function tensor and r(r, t) is 
an effective correlation time of velocity increments. Cru- 
cially, r(r, t) is the minimum of the intrinsic turnover 
time T e ddy {t) at separation r, the overall evolution time 
t, and the sweeping time r/vo. Although this result con- 
firms the sweeping decorrelation effect, we argue that the 
pair-dispersion law for zero mean-velocity ensembles at 
high Reynolds numbers is different from the t 9 ! 2 sug- 
gested by Thomson & Devenish Q. Instead, we argue 
that there are distinct ranges of power-laws i 2 , t 1 , t 6 
and then t 1 again at successively longer times, exactly 
as Thomson & Devenish argued for ensembles with large 
mean velocities. We carry out careful numerical Monte 
Carlo simulations with our diffusion model which ver- 
ify these behaviors in the model at very high Reynolds 
numbers. We also present Monte Carlo results for our 
diffusion model at the moderate Reynolds numbers em- 
ployed in current KS work, and reproduce then both the 
"i 3 -law" and "i 9 / 2 -law" that have been reported in KS at 
comparable Reynolds numbers. We thus argue that the 
current KS results in the literature are not yet probing 
asymptotic regimes and the true scaling in KS at very 
high Reynolds numbers will be the same as in our diffu- 
sion model. 

The detailed analytical derivation of diffusion models 
is presented in section|TT]of the paper, and predictions for 
their dispersion laws discussed in section lnll Our numer- 
ical methods are described and validated in section IIVI , 
and then used to obtain results for mean-square particle 
separations and other statistics. A concluding section IVl 
briefly discusses the results. 
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II. DERIVATION OF THE DIFFUSION MODEL 

In this section we present the derivations of our main 
analytical results. A reader who is only interested in 
physical conclusions and not the detailed justifications 
may skip to our final formula (|54[) for the pair-diffusivity 
and the following discussion. 



A. Gaussian Integration-by-Parts Identity 

We show first that the transition probability of par- 
ticle pairs in Gaussian velocity ensembles obeys an ex- 
act evolution equation, as a consequence of the well- 
known integration- by- parts identity or Donsker-Furutsu- 
Novikov relation (see section 4.1). Let u(x, t) be the 
random turbulent velocity field and let the fluid particle 
position that satisfies 



dt 



x(i) = u(x,t), x(i ) = a 



(4) 



be denoted as x u (a, to\t), or x(a, t) for short. Define the 
"fine-grained PDF" of 2-particle positions as 

2 

P2,u(x2,xi,£|a 2 ,ai,io) = S d (x n - x u (a„, t \tj). (5) 

Then the PDF of 2-particle positions is given by 
P 2 (x2,Xi,i|a 2 ,ai,io) = (P 2 , u (x 2 , xi, t|a 2 , ai, t )), (6) 

where the average is over the random velocity field u. 

Taking the time-derivative of (0 and using (j4]) it is a 
calculus exercise to show that 

2 

d t P2,n(t) =~Y1 Vx «' [(u(Xn,i) + U / (x„,i))P 2 , u (t)] , 
n=l 

(7) 

where the velocity has been decomposed into its mean 
and fluctuating part u(x, t) = u(x, t) + u'(x, t). The av- 
erage of the second term on the righthand side can be 
obtained using Gaussian integration-by-parts fl3j 



«(x,i))p 2jU (t)) 



J d d y J ds C ik (x,t;y,s) 

, f v A,u(*)\, (8) 
du k (y,s) I 



where dj (x, t; y, s) = (u^ (x, t)u'j (y, s)) . To represent the 
functional derivative we introduce the Lagrangian re- 
sponse function 



Gij(a,t;y,s) = 



6uj(y,s)' 



(9) 



so that 
6 



$uk(y,s 



P 2 ,u(i) = J2 -^„ p 2,u(i) • G jfc (a m ,t;y,s). 

m— 1 

(10) 



The result of averaging (JT)) is the drift-diffusion equation 

2 

d t P 2 (t) =-J2 V x „- [u*(x„,i)P 2 (i)] 

n=l 

2 

+ d < d xL [D ij (x n ,yi m ,t,to)P2(t)]. (11) 

n.m— 1 

with 



u*(x,i) = u(x,t) + a^Aj(x,x',i,i )| x , =J 



(12) 



the mean velocity plus a fluctuation-induced drift, and 
with the diffusivity tensor 



Dij(x n ,x m ,t,t ) = I ds I 

J t J 



d d y C ifc (x„,t;y, s) 



: (Gjk (a m , t ; y, s) |x 2 , Xi , t ; a 2 , ai , i ) (13) 



where 



(Gjk (a rn , i; y , s) |x 2 , Xi , t; a 2 , ai , t ) 

(G jk (a m ,t;y,s)P 2 , u (t)) 

P2(t) 



(14) 



is the conditional average of the response function given 
that the two particles start in locations a 2 , ai at time to 
and end up at locations x 2 ,Xi at time t. 

We now develop a more useful expression for the re- 
sponse function Q. It is straightforward to show by 
functional differentiation of the equation of motion (|4]) 
that 

dtGij = ^(x(a,t \t),t)G kj +6 lJ 6 d (y-x(a,t \t))5(t-s). 

(15) 

This equation may be solved as 

■fa t-v s) = i •My> s l < ) (5d (y- x ( a ' < o|s)) t> s>t 

' ,y ' ' 10 o.w. 



G t 

(16) 

with g(y, s\t) = Tcxp ( /* dr f^(x(a, t \r), r)) the time- 
ordered exponential matrix for the trajectory which sat- 
isfies x(a, to I s ) = y- This notation is made natural by an 
alternative derivation of (fl~6|) based on the flow composi- 
tion identity 



x(a,t |t) = x(x(a,t |s), s\t). 



(17) 



Taking the functional derivative 5/6uj(y, s) of (JT7J) and 
using the chain rule gives 



Sxi{a,t) dxi 
duj(y,s) dy k 



Sx k (a, s) 
y=x(a. s ) Su 3 (y,s)' 



(18) 



On the other hand, it is readily seen that the functional 
derivative of the integral form of the particle equation of 
motion (@|, gives 



Sx k (a,s) 
Suj(y,s) 



d jk S d {y-x(a,t o \s))0{s-t o ). 



(19) 
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Thus, eg.([T6|) is rcderived with gij{y,s\t) = ^(y,s\t). 
If ([To]) is substituted into the formula (|T3|) it yields 



is sizable dependence on ai , a 2 . With this motivation, we 
make the Markovian approximation 



A,(x„,x m ,Mo) ee fds [ d d y rn C ik (x n ,t;y m ,s) A;(x*, W,*o) = j ds J d d y m C ik (x n ,t;y. 



x(9jk(y m ,s\t)\x 2 ,xi,t;y m , s;a 2 ,a 1 ,t ) 
xP(y m ,s|x2,xi,t;a 2 ,ai,to) (20) 

where P(y m , s|x 2 , Xi, i; a 2 , ai, ig) is the conditional prob- 
ability density of the position of particle m at time s given 
the positions of both particles at times t and to- This for- 
mula for the diffusivity when substituted into (|TT| . (|T2|) 
gives the final form of our exact evolution equation for 
the 2-particle transition probability. 



B. Markovian Approximation 

Despite appearances, the evolution in the exact equa- 
tion dnji is non-Markovian in general. It is clear from 
formula (|20]l that the 2-particle diffusion matrix is a func- 
tion not only of the particle positions Xi, x 2 at time t, but 
also of the positions ai,a 2 at time to- This dependence 
was suppressed in our notations, but the evolution, in 
principle, retains a long-time memory of the initial con- 
ditions. Only in special cases can the evolution be shown 
to be Markovian. The famous example is the Gaussian 
velocity field that is delta-correlated in time, the so-called 
Kraichnan model flil fl5| , for which 

C lk (x, t; y, s) = C lfc (x, y; t)S(t - s). (21) 

Substituting into (|20[) and using 



9jk(y m ,t\t) = 8 jk (22) 



and 



P(y m ,i|x2,xi,i;a 2 ,ai,io) = S d (y m - x m ) (23) 

gives (with the delta- function rule" for the upper limit 
of integration) 



Dij (Xnj X m , to) — Cij[X ri 



,t). 



Thus, in this case rigorously there is no dependence of 
the diffusion matrix D on ai 1 _a 2 and the well-known dif- 
fusion model is obtained [III, Gil]- Another example with 
Markovian particle evolution is the velocity field obtained 
as the superposition of Gaussian random wave trains with 
very high frequencies, so that the group velocity of the 
waves greatly exceeds the root-mean-square velocity [l6| . 
This example has direct relevance to KS simulations with 
"eddy-turnover frequency" uj n = \y , k^ l E(k n ) in the limit 
A 3> 1 of large "unsteadiness" parameter. 

The description as a diffusion should generally hold 
reasonably well if the correlation time of the Gaussian 
velocity field is short enough, since the integrand in ([2H)l 
then becomes negligible at values of s < t for which there 



x (9jk(y m , s\t)\x 2 , xi, t; y m , s)P(y m , s|x 2 , xi, t).(24) 

The physical assumption is that for times ordered as 
to <C s < t the position of the particle at time s is deter- 
mined mainly by its position at time t and is negligibly 
dependent on the position at the initial time to- The 
worst case for this approximation is clearly the "frozen 
velocity" model with infinite correlation time, when times 
s > to in the integral are not suppressed by decay of cor- 
relations. Such s values give an undamped contribution 
also in general for times t — to much smaller than the ve- 
locity correlation time. However, it is easy to check that 
the exact result (fT3")) [or (|20[) ] and the Markovianized re- 
sult ([24]) both give 



dt 



Dij(x n ,x m ,t,ta) = Cij(x n ,t;x m ,t) + 0(t - to) (25) 



so that, for t — to much smaller than the correlation time, 

D ij (x n ,x m ,t,t ) = Cy(x„,io;x m ,to)(i-*o)+0((t-£o) 2 )- 

. (26) 

Thus the Markovian approximation becomes exact in this 
limit. We note in passing that the Kraichnan-Lundgrcn 
theory of 2-particle dispersion [l?], (3 when applied to 
the Gaussian velocity ensemble gives a result almost iden- 
tical to the formula (|24|) (for more discussion, see [lj|). 

The formula (|24p from the Markovian approximation 
can be further simplified. It is intuitively clear that con- 
ditioning on the location of both particles is superfluous 
in an average of a random variable that involves only one 
of these particles. In fact, it can be easily established 
from the definitions (O,© that 

p( y ,.s|x',x,t) - J dVP(y',y,*|x',x,t) 

= I d d y' (6 d (y> - x(x', t\s))S d (y - x(x, t\s))) 
= (<5 d (y-x(x,t| S ))) =P(y,s\x,t). (27) 
A similar argument gives 

{9jk(y,s\t)\x',x,t;y,s) = (g jk (y,s\t)\x,t;y,s). (28) 
More generally, we may define the PDF 
■P(gj t; y 1 , y, six', x, t) = 

(S dx '"U " g(y, s\t))S d (y' - x(x>,t\s))5 d (y - x(x, t\s))) 

(29) 

and mimic the previous argument to show that 

P(g, t; y, s\x', x, t) = P(g, t- y, s|x, t). (30) 

Then 

P(g,t;y,s|x',x,t) 



P(s,t\y,s;x',x,t) 



P(y,s\x',x,t) 
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P(g,f;y,s|x,f) 

= P(g,t|y,s;x,i). (31) 
It follows from these facts that 

Ai(x ra ,x m ,Mo) = J ds J d d y m C ik (xL n ,t;y m ,s) 

x (9jk(y m , s\t)\-x m , t; y m ,s)P(y m , s|x TO , i), (32) 

which is the final form of the Markovian approximation 
for the diffusion tensor. 

We now consider the special case when the velocity 
field is statistically homogeneous in space. In that case, 
the drift velocity in (|12[) is independent of x and sim- 
plifies to u*(t) = u(t), due to homogeneity and incom- 
pressibility [31] . Furthermore, a simplified equation can 
be derived for the transition probability of the 2-particle 
separation vector r = X2 — xi , defined by 

P 2 (r,t\r ,t ) = J d d a P 2 (x + r, x, i|a + r , a, i ) 

= J d d a P 2 (x,x-r, t|a + r , a, t ), 

(33) 

which is also independent of x. Since the diffusion tensor 
Dij(x n , x m , t) depends only on the difference x„ — x m in 
the homogeneous case, the equation (fill with the sub- 
stitutions r = X2 — xi and 

V X2 — > V r , V Xl — > -V r , (34) 

yields the diffusion equation 

d t P 2 (r,t\r ,t ) = d r id rj [K ij (r,t,t )P 2 (r,t\r ,t )] , 

(35) 

with the eddy-diffusivity tensor 
Kij(r,t,t ) 

= 2D i:j (0, 0, t, t ) ~ D i:) (r, 0, t, t Q ) - D ti (0, r, t, t ) 
= J ds J d d yS lk (r;y 7 t,s) 

x(g jk (y,s\t)\0,t;y,s)P(y,s\0,t) 

(36) 

and we define the 2nd-order structure function of velocity 
increments at two points 0, y and two times t, s : 

S ik (r;y,t,s) = (K(r,t)-^(0,i)][u' fe (y+r, s)-u' k (y, a)]}. 

(37) 

If furthermore the velocity field is assumed to be sta- 
tistically stationary in time, then we can take t — to — > t 
and to — > 0, to obtain 

d t P 2 (r,t\r o ,0) =8 r ,d r] [^(r,t)P 2 (r,i|r ,0)], (38) 

with 

Kij (r, t) = f dr J d d y S ik (r ; y , 0, r) 

x (g jk (y, 0; y, r)P(y, r|0, 0) (39) 

by the change of variables r = s — t. 



C. Structure Function and One-Particle 
Distribution Function 

The integral over y in the above formula (|39]l con- 
verges at large y because of decay in the two-point struc- 
ture function and in the 1-particle transition probability. 
Physically, rapid decay is due to the facts that increments 
separated by great distances are uncorrelated and parti- 
cles have low probability to be swept to large distances. 
Both of these effects can be easily quantified. 

To evaluate the two-point structure function, we use a 
standard identity that expresses it in terms of the single- 
point 2nd-order structure function ([20], p. 102): 

5 4fc (r; y, 0,r) = 1 [S lk (y + r, 0, r) + S ik (y - r, 0, r) 

-2S ik (y,0,r)}. (40) 

We first consider the single-time case with r = 0. For 
the spatial power- law covariance ([3]) with < a < 2, the 
single-point structure function becomes 

Sii(T)= 2[C^(0,r)-Cy(r ! r)]| T=0 
- 2D ir a [(d+a- 1)<% - afifj] + 0{r 2 /L 2 ){Al) 

for r <C L. The formula (|40[) implies in general that 
Sik(r\ y, 0,t = 0) ~ Si k (r) for y <C r, whereas in the 
particular case (|4"Tj) it gives 

S lk (r;y,0,r = 0) = O(r 2 /y 2 - a ) (42) 

for r -c y -C ^- When y ^> L, there is generally expo- 
nential or fast power-law decay, depending on the precise 
assumptions about the fall-off of the spectrum at low k. 
The 2-time structure function (r; 0, r) shows a simi- 
lar behavior as the single-time structure function, except 
that there is a new length Lp(r) = (D^It]) 1 ^ with eddies 
smaller than this scale decorrelated by time |t|. As seen 
from (J5J), the decorrelation is associated to an exponential 
decay of the cospectrum, with Lp(t) acting as an effec- 
tive "dissipation scale." Thus, Sjj(r;0,r) scales cx r 2 for 
r <C Lp(r), while formula (J4TJ) holds for Lp(r) «r«L. 
Thus, the decay law ([42]) is found when r^0 only for the 
range of values max{r, Lp(r)} <C y <C L and is limited 
to times |r| < L^/D^. For y <C max{r, L / 3(r)} instead 
Si k (r; y, 0, r) is independent of y and for y ^> L the de- 
cay is again like that for r = 0. 

The 1-particle transition probability should be domi- 
nated by large-scale sweeping and thus have the form 

P(y ' T| °' 0) = (27r)'/'„g|7f CXPHy " nT|2/2Uo2r2) 

(43) 

to a good approximation, with vq the root mean square 
velocity. We hereafter consider mainly the case u = 0. 
For the Gaussian random field with mean zero and co- 
variance ([2]), vo oc D 2 L a . In that case, it has been veri- 
fied by a formal scaling analysis in Q , section 7, that the 
leading-order motion of particles for large L is indeed bal- 
listic with a constant, random velocity v = y/r chosen 
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from a Gaussian ensemble with rms value Vq ■ Subleading 
corrections were also obtained in Q to account for the 
effects of the change of the velocity in space and time. 
Note that l|4*3")) decays rapidly for y 3> Vq\t\. 



rj/v and 



A(t) 



1 — a 
2(2 - a) 



1 

ja-2 



for \t\ < 1 
for \t\ > 1 



(46) 



D. Stability Matrix 

The most difficult term to evaluate in ([55)1 is the condi- 
tional average of the stability matrix g(y, r|0). Existence 
of this matrix requires a short-distance cutoff r\ on the 
"inertial-range" scaling behavior in the model covariance 
@ and ^ , which otherwise corresponds to velocity fields 
only Holder continuous and non-diffcrcntiable in space. 
Even with the cutoff, the matrix g(y, r|0) will grow ex- 
ponentially in |t| almost surely, with rate determined by 
the leading Lyapunov exponent A cx (Di/?? 2- ") 1 / 2 . It is 
thus far from clear a priori that the conditional average 
even remains finite in the limit rj — > 0. 

We begin by evaluating this term for the "frozen" ve- 
locity field with infinite correlation time (or D3 = in 
eq.©). A key observation here is that the Gaussian tran- 
sition probability (|43[) implies that particles are swept 
from point y to in time |t| along straight lines with 
a constant speed v = y/r generally of order vq. The 
velocity-gradient field Vu(x, t) has a spatial correlation 
of order 77, so that the particle trajectories contributing 
in (I39p will see a constant in space but rapidly changing 
velocity-gradient with a correlation time ~ tj/vq. Thus, 
one can expect that the Lagrangian velocity-gradient will 
be well approximated by the model of a Gaussian field 
that is delta-correlated in time, for which the statistics 
of the stability matrix has been much studied. 

To make this argument more formally, consider 
the spatial covariance of the velocity-gradient in the 
frozen case Cy, mn (r) = {u' l m (v)u' j n {0)), where u' im = 
du'i/dxm. By twice differentiating ([3]) and then averag- 
ing over the direction of the unit vector f , it is calculated 
to be 



Ci 



with D[ 



D[r a - 2 [{d- 



-(6 



Sj n + S in 5 jm )](4A) 



J 1 q(q-2) 



a— 4 
d+2 



> for d > 2 and 



< a < 2. This covariance holds for r > 77, whereas the 
covariance for r < r\ is essentially constant and can be 
taken to be given by (|4"4")l with r = 77. A particle swept 
with velocity v will see a random velocity-gradient with 
temporal correlation obtained by substituting r = vt in 
(|44| . Thus, the (Eulerian) velocity-gradients in a La- 
grangian frame can be taken as Gaussian with covariance 



(< m (tK B (0)) = Df- S n (t) [(d+ 



-(<5 



Sjn + SinSjm)] (45) 



with D" = 2 



2- a 
1-a 



D[ and 5 v (t) = f A(f ), for t n 



Since A(t) is intcgrable for a < 1 with f_ °° dt A(t) = 1, 
one then has lim^o 6 v (t) = S(t). It follows from these 
arguments that the velocity-gradient experienced by the 
particle should be approximated by a Gaussian matrix- 
valued process, constant in space and delta-correlated in 
time, if a < 1. This approximation could break down 
for fixed 77 <§C i if there happens to be a small advection 
speed 11 < Do- 
Now consider the non-frozen velocity field, with covari- 
ance given by ([2]) and ([3]) for D 3 ^ 0. In this case the 
single-point, 2-time covariance of the velocity-gradient 
averaged over directions has the form 

e«,m»(r = 0,t) = Dirf-V^MA'" [(d + l)%<5 m „ 

\^im^jn ~\~&in^jm)\ (^^) 

There is now a short correlation time t v = rf /-D3, which 
allows us to write 



(<™(*K„(o)> = -fV +/, -X(*) t( d + W 



-{SimSjn + SinSjm)} (48) 



with <5 ?; (£) = f A(J-) for A(t) = 2exp(-|£|). Thus, the 
single-point statistics of the velocity-gradient becomes 
temporally delta-correlated for vanishing 77. In addition, 
there is the same decorrelation effect of rapid sweeping 
through space that occurs in the frozen-field case. The 
latter will dominate when r\jv -C rf /D 3 at small 77 and 
when the spatial decay of correlations is fast enough, that 
is, when both /3 < 1 and a < 1. In any case, we obtain 
again a model for Lagrangian velocity-gradients that are 
Gaussian, constant in space and delta-correlated in time. 
There is here no problem with small speeds v <C vq, since 
the correlation time will never be larger than rf jD-j,. 

The stability matrix has been well-studied for Gaus- 
sian velocity-gradient fields, constant in space and white- 
noise in time. In particular, it has been shown in (2l| 
that the matrix random process g(y,r|0) is a diffusion 
on the group SL(d) of d x d matrices with determinant 
1 . We shall use specifically the formula for the transition 
probability density p T (g) of this process starting at the 
identity, Eq.(7.14) in 0] for n = 2 : 



SL(d) 



Mg)/(gr )dM(g) = / P 2 (r,r|r ,0)/(r)d d 7 



(49) 



where fi is Haar measure on SL(d) and 

9 T P 2 (r,T|r ,0) = A^ 2 P 2 (r, r|r , 0) (50) 

for 

M 2 f(r) = D[(d+l)6 l3 r 2 -2r irj }d rt d r J(r) 
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Dd ri d rj { [(d + l)5ijr 2 — 2riVj]f(r) } ,(51) where the angle-averaged structure function is defined by 



where the second line follows by incompressibility. This 
implies also that the operator is self-adjoint. Since 
M.ij = for a general linear function /(r) = a-r, and 
considering in (|49]) arbitrary choices of a, To, it follows 
that 



SL(d) 



gPr(g) d^l(g) 



(52) 



the identity matrix. This result is due essentially to the 
fact that a diffusion leaves invariant a linear profile. Fi- 
nally, we can conclude that 



(<7 ife (y,r|0)|0,0;y,r) = 5 jk . 



(53) 



The exponential growth of the individual realizations is 
offset by their rapid rotation in space which leads to large 
cancellations in the ensemble average. Incompressibility 
was necessary to the argument. 

The result ([5U]) is only strictly known to be valid when 
the velocity covariance converges to an 77- independent re- 
sult as 77 — > 0, whereas (|45|) diverges as ~ rf 1-1 for a < 1 
and (|4"51) diverges as ~ f Qr q, _|_ ^ < 2. However, 

the final result (|53[) is independent of the amplitude of the 
covariance (i.e. the value of D\) and thus we conjecture 
that it extends even to the present cases with diverging 
covariance. The result yields a simplified formula for the 
2-particle eddy-diffusivity: 

K ij (v,t) = J dr J d d yS ij (r;y,0,r)P(y,T\O,0), (54) 



together with (|40l) . (|43|) . We shall now use this formula 
to obtain concrete results for the eddy-diffusivity in the 
Gaussian ensembles whose covariances are given by 



Sij{r;y) = / dOAi( r ;y) 



(58) 



for S d = 2ir d / 2 /T (f) the (d - l)-dimensional area of the 
unit hypersphere in d-dimcnsional space. 

When the velocity statistics are isotropic, as for the 
model with zero mean and covariance ([5]), the eddy- 
diffusivity tensor can be reduced to two scalar functions 
Kl,Kn defined by 

Kij{r, t) = K L (r, tjnfj + K N (r, t)(% - f t fj). (59) 

These two functions are related by incompressibility as 
A'tv = Kl + rK' L /(d — 1) and it is convenient to base 
further analysis on Kl ■ As is well known, if the separation 
statistics are also isotropic, then the diffusion equation 
([331) can be expressed entirely in terms of Kl , as 



1 d 



dP 



d t P ^) = —Q- r [r d - l K L {r,t) — {r,t) 



dr 



Here the separation PDF satisfies 

P{r,t) r d - x dr = 1. 



(60) 



(61) 



as normalization condition. 

The displacement vector y in (|55j) breaks rotation 
invariance, but the average over solid angle restores 
isotropy. We can thus decompose also 

Siji^v) = SL{r;y)rirj + S N (r;y)(Sij -hrj) (62) 

into longitudinal and transverse contributions with re- 
spect to the separation vector r. By dimensional analysis 
one can write 



E. The Frozen-in-Time Velocity Field 



S L (r;y) = S L (r)F 



S L (r)F(^j, (63) 



The simplest case to analyze is the "frozen" field so 
that 

S,(r;y,0,r)=S y (r;y). (55) 
Making the change of variables u = y 2 /2v 2 t 2 , 



drP(y,r|0,0) 



1 fd-l y 2 

-r ' 



V8ir d / 2 voy^ 1 V 2 '2v$t 2 



(56) 

with the (upper) incomplete gamma function defined by 
T(s,z) = l™duu s - x e- u . Since d d y = y^dydfly, with 
dQ y the element of d-dimensional solid angle, we get from 
(1541) that 



K a ( r > *) = /or 7T\ / — Si i ( r ' y) T 



V2Y (i) 7 "0 



2 ' 2v 2 t 2 

(57) 



the latter for L ^ r. The function F{y/r) can be in- 
terpreted as the correlation coefficient of (longitudinal) 
velocity increments 5vl(t) at points a distance y apart. 
For the velocity covariance $2$ with D 3 = it is possible 
to derive a complicated, closed-form expression for the 
function F(w) as suitable combinations of Gaussian hy- 
pergeometric functions of the argument w 2 . However, we 
shall not pursue this here. The most important property 
of F, which follows from (|4"2")l , is 

FH ~ I ) ^ (2 a ) W < ! ■ (64) 
v ; \ {const. )w- {2 ~ a > w > 1 v ' 

Thus, we can write KL(r,t) = Sl{t)t{t, t) where 



r(r,t) 



V2r(|)7 vo Vr 



dy_ F fy\ T [d- 1 y 



2 '2v 2 t\ 



(65) 
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is a 2-particle Lagrangian correlation time. With the 
substitution y = rw, this becomes 



Since the y-integration has the form of a convolution, it 
is easily evaluated by a Fourier transform: 



r(r,t) = — J(x), x = 

v r 



(66) 



for 



JW = 7iF(I)f ^K^'^) (67) 

For d = 3, r(l, z) = e~ z and (|67|) is a Laplace transform 
in the variable w 2 . The most directly useful consequence 
of ([57)1 is the asymptotic behaviors 



J(x) 



X I«l 
Joo X > 1 



where we have used J °° dv T 



we have defined 



d-i 



= } 2 l I dwF(w) 



(68) 



\/2Wf) and 



(69) 



The latter integral converges for a < 1 . We conclude that 



r(r,t) 



(70) 



Our result is quite similar to that obtained by |8[ for 
the case of large mean velocity u; see their equation (8). 
Some differences are that our eddy-diffusivity is isotropic 
and has the short-time behavior proportional to t. How- 
ever, most importantly we see the same sweeping decor- 
relation effect, with the 2-particle correlation time at long 
times proportional to the sweeping time t/vq. With no 
such effect one would instead expect the correlation time 
to be always proportional to t in the frozen-field case. It 
should be emphasized that we obtain this result in the 
zero mean-velocity ensemble, where Q have predicted 
different behavior. We shall compare our results with 
theirs in more detail in section IIII1 where we shall also 
derive the quantitative predictions of our formula for the 
growth of mean-square particle separations. 



F. Finite Time-Correlated Velocity Field 

We now study the Gaussian model with covariancc 
([2]) for L>3 7^ 0. More generally, consider any velocity 
field statistically homogeneous in space and stationary 
in time. Then (J32]) together with (gSJ) & (J53J) give 



Dij (x' - x, f ) = J dsj d d y Cy (x' -y,t-s) 

|y-x-u(i- S )| 2 l 



[2lTV 2 (t - s) 2 ] d /' 



■ cxp 



2vl\t-s\* 



(71) 



Dij(k,t) = / dsdjfct-a) 
'o 

1 



x exp 



ik-H(t - s) - ^vlk 2 (t- s) 2 



(72) 



For the model in @ note C7y(k, t) = Cij (k) exp(— 7fc|i|) 
with C, u (k) = D 2 P ij (\s)/k d L +a and 7fe = D 3 k P L . For large 



u, see [S|J. Hereafter we take u = 0. Then making the 
change of variables a = vok(t — s), one obtains 



1 



v k Jo 



Vokt 

da Cij (k) 
x exp 



\v kj 



1 



Thus, for t < 1/vok, 

A 3 '(M)~Cy(k)t. 
This implies by inverse Fourier transform that 
ify(r > t)~S , y(r)t > t<r/v . 



(73) 



(74) 



(75) 



On the other hand, consider fixed t and large k. Note 
that the convection time is smaller than the correlation 
time, or vok > 7^, for k > fc* = (D^/vq) 1 ^ when j3 < 1. 
Thus, for k > fc*, (73]) gives 



D«(k,t)~— C^(k) M /-erf 



75 



(76) 



This formula is exact in the case of frozen turbulence 
(-D3 = 0) when fc* = 0. If furthermore k 1/vQt, then 



(77) 



becomes independent of t and scales as a power 
^,-(d+o:+i)_ p or a < wc thus obtain by inverse Fourier 
transform that for r <C minjuot, 



7rD (a+1) 



(78) 

It follows that the essential behavior of the frozen field 
case carries over to the finite time-correlated velocity 
with a < 1 and /3 < 1. Just as for the frozen velocity, 
Kl(t, t) = <Sx(r)r(r, £) and the correlation time satisfies 

([SB] and ASS]) with = J\ 2^-[H. 

If, however, /3 > 1, then the behavior is quite different. 
Under this assumption 7^ > «ofc for k > fc», so that we 
now make the change of variables a — 7^(i — s) to obtain 



A,(k,i) = i 



7fc Jo 



7fc* 



dtrCy(k) 
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x exp 



1 f vok 



2 V Ik 



(79) 



Equations (|74|) and (|75|) again hold, now for t l/"fk 
and t <C l/71/rj respectively. On the other hand, for 
fixed £ and fc 3> A;* = !/£*, 



Dy(k,i) Cy(k) [1 - exp(- 7fc t)] 



7fr 



If furthermore k S> 1/Lp(t), then 

A;(k,t) ~ — Cy(k) 
7fe 



(80) 



(81) 



becomes independent of t and scales as a power 
fc-id+a+p) _ when q- + /3 < 2, we then obtain by inverse 
Fourier transform that for r <C min{L ( g(t), L*} 



(«+/3) 



-r oW [(d + a + /3)% - (a + + 1)^] 

(82) 



We can again write Kl(t, t) = SL(T)r(r,t) but now 



r(r,t) 



t t < r^/Da 

(const. t > r^/D 3 , r < X* 



(83) 



Thus, the sweeping decorrelation effect is absent at suf- 
ficiently small scales when (3 > 1 and a + /3 < 2. 



III. CONSEQUENCES OF DIFFUSION MODEL 

In the previous section we have derived a diffusion 
model which, for homogeneous and isotropic statistics, 
takes the form ([50]) . For the Gaussian velocity ensemble 
having covariancc @ with Kolmogorov scaling exponent 
a = 2/3, the diffusivity takes the form 



K L (r,t) 



CL£ 2/3 r 5/3 



J 



Vo 



r 



C L (er) 2 / 3 t t-^r/vo 
, 2/3 5/3 



(84) 



both in the frozen case and in the temporally fluctu- 
ating case with j3 = 2/3. Here Cl is the Kolmogorov 
constant in the longitudinal velocity structure function, 
S L (r) ~ C L (er) 2 / 3 , and C' L = ClJ^. In this section 
we shall attempt to determine the growth law for the 
mean-square separation (r 2 (t)) predicted by the model 

Does this model lead to the £ 9//2 -law of Thomson- 
Devenish To answer this question, we must briefly 
review the argument for the 9/2-law. The key idea 
in Q is that the mean-square separation pointwise in 
space depends on the local value v' of the fluctuating 
velocity. The sweeping effect occurs at points where 
T S w(f) = r / v ' is smaller than the intrinsic correlation 



time, Tint(r,t) = e _1 / 3 r 2 / 3 for finite-correlated velocity 
(f3 = 2/3) and r int (r,t) = t for "frozen" velocity. The lo- 
cal correlation time is argued to be the smallest of these: 



-{r,t) = mm{T aw (r),Ti n t(r,t)}. 



(85) 



Hence, when v' > (er) 1 / 3 (fluctuating) or r/t (frozen), 
then the mean-square separation conditioned on v' is af- 
fected by sweeping and shows the slow growth 



'(*)>* - — 



(86) 



but in the opposite case exhibits the faster growth 

(r 2 (t)) v , ~et 3 . (87) 

Using these growth laws to evaluate t sw and r int in l[55|). 
it is easily checked that the i 6 -law holds for points with 
v' > (et) 1 / 2 and the i 3 -law for points with v> < (et) 1 / 2 . 
The probability for the latter condition to hold is small 
but growing in time: 



Probfi/ < (et) 



1/2 



(et) 



3/2 



This formula holds for a Gaussian distribution of 3D 
velocities v', or for any similar distribution p(v') = 
(l/i>g)/(v'/z;o) with variance v 2 , and non-vanishing den- 
sity at the origin. The unconditional mean-square sepa- 
ration can then be estimated from ([57)) and (|55|) as 



et 3 



(et) 



3/2 



5/2^9/2 



(89) 



The same result can be obtained from the t 6 dispersion 
law (|86|) by noting that it is a rapidly decreasing func- 
tion of v' , so that the dominant contribution is obtained 
from the points with v 1 > (et) 1 / 2 which also occur with 
probability ~ (et) 3 / 2 /vq. 

At first sight, it appears that the model (|60 p .(|84 p may 
embody these ideas of §. The diffusion model implies 
the exact equation 



l(r 2 (t))=2 J K T (r,t)P(r,t)r d -Hr, 



(90) 



where Kt = Kl + (d — 1)Kn is the trace of the diffusion 
tensor. The average over r in (f9"0")) can thus play the same 
role as did the average over v' in the argument of Q . The 
eddy-diffusivity ((54")) is equivalent to the correlation time 
(j?0")) . The population of particle pairs with separations 
r > i>ot should exhibit a growth law (?" 2 (i))> ~ et 3 , while 
the pairs with r < v^t should exhibit (r 2 (t))< ~ It 

v 

appears possible that averaging over the entire range of 
pair separations could give rise to the 9/2-law ([55)) with 
an intermediate growth rate. 

The above reasoning is, however, essentially wrong. 
The diffusion model ((6T))) . (|84[) does possess a t 3 regime, 
but only in an unphysical way. To see this, note that for 
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both the t 3 and the t 6 growth laws the condition r > vot 
is first satisfied only at such long times that t > Vq/s. 
Substituting the standard relation e ~ Vq/L (which fol- 
lows from the assumed Kolmogorov scaling of the energy 
spectrum) implies that the t 3 law can be self-consisently 
satisfied only for times greater than a large-eddy turnover 
time, t > L/vq. In that case, (r 2 (t)) exceeds L 2 and the 
particle pairs have left the inertial range. As we shall ver- 
ify below, the model (|60|) . (|84|) docs indeed possess a t 3 
range when t > v 2 /e but this exceeds the validity of the 
model, which was derived only for the range r < L. In the 
range r > L the two particles should instead execute in- 
dependent Brownian motions with a constant diffusivity 
Dt and the mean-square separation grow diffusively as 
(r 2 ) ~ ADrt. Thus, the t 3 range is an unphysical artefact 
of the model ([50 ]) . (jgl j) . 

The argument for the asymptotic t 9 / 2 law by Thom- 
son & Devenish thus fails for the model ([60)) . (|84|) . 
Very importantly, however, we shall show below that our 
diffusion model can produce an apparent t 9 / 2 law over 
a finite range of scales at relatively low Reynolds num- 
bers, for similar choices of parameters with which such 
growth laws have been observed in kinematic simulations 
[3, [HI [ll| . In this case, the eddy-diffusivity in the equa- 
tion (|60[) is not given by the formula (|84[) , which is asymp- 
totically valid only for L»r, but instead directly from 
the expression ([75| . which holds in general. We shall thus 
suggest that the 9/2 growth law observed in several kine- 
matic simulations is a finite-Reynolds-number effect and 
does not represent the asymptotic behavior that would 
be observed with very long inertial ranges. 

We argue that the true high-Reynolds-numbcr 
behavior — both in our diffusion model and in the kine- 
matic simulations — is essentially the same as that found 
by Thomson & Devenish for the situation of large mean 
velocity u section 3.1). The principal difference 

is that we obtain also an early-time Batchelor ballistic 
range (22L [23j with t 2 growth. This is followed, as ar- 
gued in [8j, by ranges of diffusive t 1 growth, t e growth 
and finally by a range of t 1 or t 3 growth, depending upon 
whether the correct diffusivity (fT3")) is used for that range 
or whether the r <C L approximation (|84[) is used (in- 
appropriately, since r ^> L). We have not been able to 
find an analytical solution of our model ([BUI) . ([8"4"]) which 
exhibits all of the above ranges. In this section we shall 
instead argue using a simple mean-field approximation 



dt 



r 2 = 2K T (r,t), r(0) 



(91) 



which ignores fluctuations in the random separation r. 
In the following section [IV] we shall verify our theoretical 
conclusions by a numerical Monte Carlo solution of the 
diffusion model. 

The Batchelor t 2 regime is the only one which we can 
derive directly from our diffusion model (|60p without any 
approximation. We take as our initial condition for the 



diffusion equation the spherical delta function 

S(r - r ) 



Po{r) 



„d-l 



(92) 



with all pairs initially at separation r$. If this is substi- 
tuted into the exact equation ([90]) , it yields 



s<' J < f » 



0, 



2SrM (93) 



t=o 



where the trace of the short-time result (|75p was used, 

K T (r,t) ~ 2S T (r)t, (94) 

for t <C t~o /vo • Taylor series expansion then gives 

(r 2 (t))-r 2 = S T (r )t 2 +O(t 3 ), (95) 

which is the well-known result of Batchelor 0, [23| . The 
mean-field approximation (|91l) is exact in this regime, 
since sufficient time has not passed to change r substan- 
tially from its initial (deterministic) value r . 

As noted in Q , there is an interval of times t > r^/vo 
when r has still not changed substantially from its initial 
value r . For r « r but t ^> r /v a , the result (|M|) is 
replaced with 



K T (r 0l t) - K T (r , 00) 



e 2 / 3 rn /3 
C' T (96) 



where C' T = ^C' L as a consequence of incompressibility. 
The growth law then becomes diffusive 



(r 2 (t))-r 2 ~2X T (r ,ooK 



(97) 



this period lasting until the "takeoff time" t to when 
K T (r , oo)t to ~ rl, or 



1/3 
.2/3 ' 



(98) 



See [a]. Together with the previous Batchelor regime, this 
diffusive range is obtained from the mean-field model (|91[) 
simplified to dr 2 /dt — 2KT(ro,t). It is interesting that 
the diffusive behavior (|9T|) at early times is the an alog ue 
in the Kraichnan white-noise advection model [3 Il5l | of 
the Batchelor ballistic range (e.g. see (24[, section II. B). 
This is not an accident. The large-scale sweeping of par- 
ticle pairs through stationary eddies produces an effec- 
tive small correlation time ro / i>o which makes the veloc- 
ity field appear to be temporally white-noise for times 
t tq/vq. This is closely connected with previous at- 
tempts to simulate the Kraichnan white-noise ensemble 
by sweeping fixed large-scale velocity fields rapidly across 
the computational domain (25l . |26| . 

For times greater than the "takeoff time" t to but 
smaller than the "end-of-sweeping time" t es = v 2 /e, one 
must solve the mean-field equation (|91[) with 



dr 2 /dt = 2C' T 



-2/3 5/3 



(99) 
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which leads to the t 6 -\aw (p}. Instead for t > t es one 
must solve 



dr 2 /dt = 2C T (er) 2/3 t, C T = ^-C L 



(100) 



at least for the model As previously discussed, this 
leads to the Richardson i 3 -law but in an unphysical way, 
since r > L lies outside the validity of the model ([84)). 
For t > t es and r > L in reality Kx(r,t) ~ 2D?, where 
Dt is the 1-particle diffusivity of Taylor [27T ]. Thus, one 
must solve 



<2r 2 /* = 4L>t 



(101) 



which yields the very long-time diffusive range. 

Our picture of particle dispersion in the zero-mean syn- 
thetic turbulence ensembles is thus very close to that in 
the large mean-velocity ensembles. This is in contrast to 
Thomson & Devenish [8j , who argue for a distinct behav- 
ior of particle dispersion in the two cases. To understand 
better why we reach a different conclusion, it is useful to 
redcrive our results for the eddy-diffusivity in a slightly 
different way. For convenience we consider only the case 
of frozen velocity fields. Taking the longitudinal compo- 
nent of the formula (|54|) yields 



K L (r,t) = j dr j d d yS L (r;y)P(y,T\O,0). (102) 

As discussed in section Hi CI the factor P(y,r|O,0) arises 
as the density of the Gaussian large-scale velocity v = 
y/r. Changing to this variable in the above integral yields 



K L (r,t) = / K L (r,t\v)exp 



d d v 



2vlJ (2jTV 2 ) d / 2 



with 

and 

r(r, t\v) 



K L (r,t\v) = S L (r)r(r,t\v) 



drF 



v\t\ 



t i < r/u 
loo 1 ; t^r/v 



(103) 



(104) 



(105) 



for loo — Jq dwF(w). Since F(v\r\/r) is the correla- 
tion coefficient of increments Su(r) at distance v\t\ apart, 
Kl{v, t\v) and r(r, t\v) can be interpreted as pair diffu- 
sivity and correlation time for given large-scale velocity 
magnitude v. It is easy to average these quantities over v 
and recover the previous results for K^r^t) and r(r, i), 
in particular formula ([701) . an d our predictions in this 
section for (r 2 (t)). It was already observed in Q (sec- 
tion 3.2, p. 292) that averaging the pair-diffusivity over 
the large-scale sweeping velocity would lead to the i 6 -law 
also for the zero-mean velocity ensembles. Thomson & 
Devenish argued, however, that correct results should be 
obtained by averaging (r 2 (t)) v rather than by averaging 
Kl{t, t\v). We find that the opposite is true. The exact 
integration-by-parts identity for Gaussian velocity fields 
leads to our formula (|103|) in which the effective diffusiv- 
ity is indeed averaged over large-scale sweeping velocity. 



IV. NUMERICAL SIMULATIONS 

We now present numerical results for the diffusion 
models derived in the previous sections, both to confirm 
our theoretical predictions of their behavior and to obtain 
new conclusions where no analytical results are available. 



A. Methods and Tests 

As in [H, we shall solve the diffusion equation ([38]) 
using a Monte Carlo method for the equivalent (Ito) 
stochastic differential equation 



dn = bij (r , t)dWj (t) , i, j = 1, 



(106) 



where Einstein summation convention is used, Wj(t) is 
a vector Wiener process and 2Kij = bikbjk, with lower- 
triangular square-root b^ calculated by Cholcsky decom- 
position. We can integrate the stochastic equations (|106[) 
using the standard Euler-Maruyama scheme: 



Ti{t k ) = n(t k -i) + bij(r,t k -i)V AtNk.j i,j = l,...,d 
t k = i fe _i+At (107) 

where N k j for j = l,..,d, k = 1,2,3,... is an indepen- 
dent, identically distributed sequence of standard normal 
random variables. The normal random variables are ob- 
tained from uniform pseudorandom numbers generated 
by the Mersenne Twister algorithm [? ] which are then 
transformed to normal by the Box-Mullcr method [? ] . 

Unfortunately, the ranges of time that we must cover 
are so large that it is completely impossible for us to 
use a constant timestep At. Instead we use an adaptive 
scheme similar to that of Q. The stepsize is determined 
over geometric intervals T(m) < t < T(m + 1) with 

T(m) = A exp(Bm) for m = 1, 2, M. (108) 

The constants A, B and M are chosen for an initial par- 
ticle separation ro as 



AI 



B 



250 



ln(10) 
ln(10 9 /r ) 



lnr + 2251 



M - 2 
A = 10^ 5 r exp(-B) 



(109) 

(110) 

(111) 
(112) 



so that T(l) < t sw , T(L) > t es and AT = T(m + 1) - 
T(m) t. In each such interval we take 



At = CArriin 



K T (r,t) 



AT 



(113) 



where Kt is the trace of Kij . A large number S of inde- 
pendent samples of the process (| 106[) are generated with 
initial separations r(t = 0) = r uniformly distributed 
over a sphere of radius |ro|, and statistics obtained by 
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averaging over realizations. Most of the results presented 
below used S = 10 4 . 

There is considerable debate in the literature, how- 
ever, whether such adaptive time-stepping schemes lead 
to converged, unbiased results for the statistics 0-|TT| ■ 
To test our numerical methods, we found it useful to 
consider somewhat simpler diffusion models where exact 
analytical results are available for comparison. The mod- 
els with a power-law diffusivity 



K L (r)=Dr c , 0< C < 2 



(114) 



have been very well studied. It has been shown that the 
long-time evolution is self-similar, with a dispersion law 



(r 2 (t)) ~ g(Dt) V \ 



7 



9 = 



4/ 7 p ^d±2^j 



(115) 



and a stretched-exponential PDF 
1 



P(r,t) 



( r 2(f))d/2 

where 7 = 2 — £, 



exp 



((r 2 (t))i/2) +/3 



r((d + 2)/ 7 )r /2 



/3 = ln 



r(d/ 7 ) 



(r(d/ 7 ))( d + 2 )/ 2 



(116) 



(117) 



(118) 



with the normalization condition J °° r d ~ x P(r,i)dr = 1. 
See [? ], eqs. (3.14), (3.22) and the general, self-similar so- 
lutions found in [12( for the case I = 33[ • Incidentally, 
note that the mean-field equation (|9Tj) leads to power- 
law growth with the same exponent 2/7 as in (|115|) but 
with a different prefactor g MF = ( 7 (rf + C)) 2 ^ 7 than g. 
It is not hard to show that g MF > g, with g MF — s- g as 
d — > 00 from Stirling's approximation. 

Notice that the inertial-range model (|8~4")l reduces to 
the time- independent diffusivity 



K L (r, 00) ^C^r^/vo, 



(119) 



as long as r <C i>oi- This is a special case of the power-law 
diffusivity (fTT4|) with C = 5/3, or 7 = 1/3, so that the 
mean-square separation grows as t 6 . This case is thus 
most suitable to test our numerical methods. For the 
purposes of comparison in the next section with the more 
complex model ([8lj). we take D = C' L e 2 ^ 3 /v with C' L = 
1.262 and d = 3 so that 



'(*)} 



F 4 f 6 
PD e 1 



(120) 



with the power-law diffusion model predicting g PD = 
15.968. This model also has the self-similar PDF of form 
(fTTrJl) with d = 3, 7 = |, so a = 11.3714, (3 = 10.1767. 
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FIG. 1: Monte Carlo results for (|r(t)| 2 ) in the power-law 
diffusion model ([119} with 5* = 10 4 and C A = 1. 



We now employ the numerical scheme discussed earlier 
to see which of these exact results we can successfully 
reproduce. As we see in Fig. Q] long ranges of perfect 
t 6 power-laws can be obtained in log-log plots. On the 
other, Fig. [5] is a semilog plot of the dispersion com- 
pensated by the analytical result (|120[) . It shows that 
the prefactor is poorly calculated by our Monte Carlo, 
which gives g§ IC = 26.10 and does not appear to con- 
verge to the analytical result g FD as Ca is decreased. 
Finally, Fig. [3] shows the logarithm of the PDF of pair 
separations r plotted versus r 1 / 3 at 14 different times 
in the long i 6 -rangc. Self-similarity is well-confirmed by 
collapse of rescaled curves for different times, but the an- 
alytical result (|116p is not very accurately reproduced. 

Our conclusion from these exercises is that the adap- 
tive time-stepping scheme should be adequate for expo- 
nents of dispersion power-laws, but not for prefactors or 
PDFs. Since the primary issue in this work is the expo- 
nents, we shall employ the adaptive schemes when neces- 
sary to cover extensive ranges where constant time-steps 
are unfeasible. As additional checks on our numerical 
results for exponents from adaptive schemes, we test for 
convergence using constants Ca ranging from 1 to 10 -6 . 
We also compare our Monte Carlo results for the diffu- 
sion equation with a separate numerical solution of the 
mean- field equation (|9Tj) , integrated with a Fortran 90 
implementation of the Watt and Shampinc RKF45 ODE 
solver [2^, [2^| . This standard ODE integration method 
is also adaptive, but with variable time-step determined 
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B. The Inertial- Range Model 
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FIG. 2: Monte Carlo results for (|r(t)| 2 ) in the power-law 
diffusion model (|119[) with S = 10 4 and various Ca , compen- 
sated by the analytical result 1(1200 . 
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3: Logarithm of the rescaled PDF of pair-separations at 
14 different times in the t 6 range, for the power-law diffusion 
model fTT9]) . Monte Carlo results for 5" = 10 5 and Ca = 
1. The straight line marked with circles (o) is the analytical 
result (fTTo) for d = 3 and 7 = 1/3. 



by preselected error tolerances. We therefore can have 
confidence that the numerical results for the mean-field 
theory are well converged. 



We consider first the model (|84|) obtained for Kol- 
mogorov scaling exponents in the limit L> 1 and thus 
physically applicable only for separations r in the inertial 
range of scales. This diffusion model applies for both the 
frozen velocity case and the finite-time correlated case 
(since (3 = 2/3 < 1). For the purpose of simplifying 
the numerical work, we opted not to use the exact scal- 
ing function J{x) given by integral (|67[) . which in three 
dimensions yields a complicated expression in terms of 
generalized hypcrgcometric functions. Instead, we built 
a function with the same asymptotic behaviors (|68[) as 
the true J(x). We took 



J(x) = Jooerf(Ax) 



x i«l 

Joo X » 1 



(121) 



with A = 



and J a 



D 



(5/3) 



= 0.6396. Our ex- 



2 D f/3) 

pectation was that only these general features should be 
sufficient to observe the scaling regimes predicted in the 
previous section. This idea was borne out by the nu- 
merical results. In Fig. 2] we plot (|r(t) — i"o| 2 ) for the 
inertial-range diffusion model with ro = 10~ 20 . On the 
same graph we plot for comparison the numerical solu- 
tion r 2 (t) — r 2 , of the mean-field equation (19"TT) . The two 
agree very well, and clearly exhibit the four predicted 
regimes with power-laws oc t 2 , , t 6 and t 3 , successively. 
A convergence analysis of our adaptive scheme for these 
results is presented in Appendix A. 

To further test the theoretical predictions, we investi- 
gate the crossover times between the different regimes 
and the prefactors of the scaling laws. For example, 
in Fig. [3(a) we show for various values of rg/L the 
quantity (\r(t) — i"o| 2 ) compensated by the Batchelor- 
range prediction -^-Ci(ero) 2 ^ 3 t 2 plotted versus the time 
t/r sw rescaled with the sweeping time t sw = tq/vq. The 
Batchelor prediction fits the Monte Carlo data to within 
0.13% relative error and the end of this regime is very 
close to t/r sw = 1. Wc similarly show in Fig. EJb) 
for the same choices of Tq/L the mean-square separa- 
tion (|r(t) — ro| 2 ) compensated by the Kraichnan-likc 
"diffusive-range" prediction ^C' L e 2 ^ 3 r^ 3 t plotted ver- 
sus t/tfo with the "takeoff time" t to given by equation 
(f9"5]) . The diffusive-range prediction is verified with a 
1.8% error and the end of this regime quite convincing 
scales as ~ 10 — 2 t to . In Fig. [SJc) we show the corre- 
sponding plot of mean-square separation compensated 
by £ 4 t 6 /i>g versus tjtto- We see that a t 6 range begins 
at time ~ 10 2 t to and extends to the end-of-sweeping 
time t es = Vq/e with a prefactor g^ IC ~ 22.96 of the t 6 - 
law. This Monte Carlo value lies between the mean-field 
prediction g£ IF = 55.61 and the exact power-law diffu- 
sion model prediction g§ D = 15.97, but is quite close 
to the Monte Carlo result g£ MC = 26.10 for the latter 
model. This suggests that the adaptive time-integration 
scheme underestimates the effects of diffusion, whereas 
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FIG. 4: Numerical results for (\r(t)— i"o| 2 } in the inertial-range 
model (|84|l : Monte Carlo solution of the diffusion equation (o) 
with Ca = 1, S = 10 4 and mean-field approximation (— ). 



the true value for the inertial-range model (|84[) is proba- 
bly <76 = 15.97, the same as the power-law model (|119p . 
It is interesting that the transition between the t 1 and 
i 6 scaling ranges is quite broad, covering about four 
decades. We show finally in Fig. [SJd) the mean-square 
separation compensated by the Richardson prediction et 3 
plotted versus t/t es . For t > t es there is a clear t 3 regime 
with Richardson constant g^ c ~ 8.97. This constant 
again lies between mean-field predictions and exact re- 
sults for a self-similar solution and should not be regarded 
as accurate. Of course, as emphasized earlier, this entire 
regime of the inertial-range diffusion model is unphysical 
and will not be observed in KS model simulations. 



C. Comparison with KS Models 

Our derivation of diffusion model approximations was 
sufficiently general that we can consider cases of more di- 
rect relevance for KS simulations, with any energy spec- 
trum and without the approximation of large L. Using 
the formula (|76[) . which is exact for frozen turbulence, 
one obtains by inverse Fourier transform in 3D that 



It is convenient to assume statistical isotropy, so that 

E(k) 




t/tto 



as 


1.15 


o 






1.1 


UJ 


1.05 


to 






1 








0.95 


U 




1 


0.9 


u 






0.85 




0.8 




1.15 








1.1 








1.05 




1 


O 




Sh 




1 


0.95 


U 






0.9 




0.85 




0.8 











(c) 









10" 



10' 



t/tto 



10' 



10" 




FIG. 5: Monte Carlo results for (\r(t) - r | 2 >, Ca = 1, S = 
10 4 . Each panel shows the same curves with different scalings. 
(a) Batchelor regime, (b) Kraichnan regime, (c) t G regime, 
(d) Richardson regime. The initial separations are ro/L = 
10- 5 (o), l(r 8 (n) , 10- u (>), 10" 14 (V), f0~ 17 (o), f0- 20 (<). 
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where Pij (k) is the projection operator onto the subspace 
orthogonal to k. The trace of the diffusivity tensor be- 
comes 

„ , „ V27T f°° dk . „ / V()kt\ smtkr) , 

D ^ t) = — l t ^ (~75 J ~^ ,124) 

and Dz,(r,t) can be recovered from 

D L (r,t) = — d l D T (p,t) p^dp. 
r Jo 

Finally, the diffusivity that appears in equation (|60j) is 

K L (r,t)=2(D L (0,t)-D L (r,t)). 

To apply these results to the KS models [l-HH , let us 
recall that those models have a discrete set of wavenum- 
bers distributed as 




(125) 



for n = 1,...,N where k\ = 2tt/L, fc;v = 2Tr/rj and 77 is 
the analogue of the Kolmogorov dissipation length. The 
energy spectrum generally adopted in these models is 

E{k) = C K s 2/3 k n 5/3 S(k - k n )Ak n (126) 

n=l 

where Ak n = (fc„ +1 — k n -i)/2 and Ck = 1-5 is the Kol- 
mogorov constant, so that Cl = 1.973 Here e is a 
constant with dimensions of energy dissipation per mass 
chosen to prescribe values of the rms velocity: 

/ 2 f^ N 

v = J - ^ E{k)dk. (127) 
The formula ([124]) with the KS spect rum (|126[) yields 



D T (r,t) = — V2n 1 sin ( k nr) Ak n 

(128) 

The assumption of isotropy in this formula is only ap- 
proximately valid for KS simulations. It would be pos- 
sible to use the general result (|122[) . without assum- 
ing isotropy, which would lead to a discrete sum over 
wavevectors rather than wavevector magnitudes. How- 
ever, this would make numerical implementation a bit 
more difficult, without essentially different physics. 

We now present simulation results for diffusion models 
based on Gaussian velocity fields with the spectra of KS 
models, or, to be brief, "KS diffusion models" . The same 
Monte Carlo method was employed as for the inertial- 
range model. In all of our numerical studies we take vq = 
L = 1. We tried various values for the number of modes 
N and we found that the numerical results on dispersion 
laws in log- log plots for N > 100 are not significantly 



different (see Appendix B). All of our presented results 
are for N = 500, a comparable number to that in the 
KS studies [sT-flll]. We have also followed the practice 
in the KS literature of choosing the smallest length-scale 
77 = ro/10, f° r initial separation r$. 

Our first set of numerical experiments investigated 
whether these more realistic models would exhibit the 
power-law scaling ranges predicted in section Hill with a 
k^/ki sufficiently large. In Fig. [6] we plot the numerical 
results for the mean-square separation (|r(i) — ro| 2 ) ob- 
tained from the KS diffusion model with ro = 10 -20 . We 
observe very clearly the predicted ranges with power-laws 
t 2 ,t x ,t® and, lastly, the diffusive t 1 range at long times 
expected for a model with finite L. For comparison, we 
also plot numerical solutions of the mean-field equation 
(f9Tjl using the diffusivity (|128|) . As before, the mean- 
field theory predictions are quite close to the full Monte 
Carlo solution of the diffusion model. Lastly, we plot 
the solution of the mean-field equation for the inertial- 
range large- L diffusivity, with the same choice of con- 
stants L, vq and e. As expected, the dispersion law from 
this approximation agrees quite well with that of the KS 
diffusion model for r < L, but predicts a spurious t 3 
power-law range for r > L. The good agreement justi- 
fies a posteriori our simplification of the scaling function 
J(x) in section IIV Bl Our most important general con- 
clusion from this set of experiments is that the KS diffu- 
sion models and, we believe, the KS models themselves 
should exhibit the above four scaling ranges with succes- 
sive power-laws t 2 , t 1 , t 6 and then t 1 again, whenever the 
scale ratio k^/ki is sufficiently large. 

In order to discriminate between various alternative 
theories, it is useful to compare predictions not only for 
mean-square separations but also for the full probability 
density P(r,t). Although we do not expect our adaptive 
time-stepping algorithm to be sufficient to reproduce ac- 
curate PDFs, it is still useful to present a few numerical 
results for the KS diffusion model. In Fig. [7] we plot 
the Monte Carlo probability distribution calculated for 
36 different times spread within the t 6 range. These are 
rescaled to test for self-similarity and collapse quite well. 
It should be emphasized that the overall evolution of our 
IR and KS diffusion models is not self-similar, globally in 
time. This can be seen most clearly in the existence of 
time ranges with distinct power-law growth laws, whereas 
a truly self-similar evolution should have just one power- 
law. In a sufficiently long t 6 range, however, one should 
expect a self-similar evolution. For example, the inertial- 
range model (|84|) in the t 6 range reduces to the time- 
independent diffusivity K^(r, 00) = C' L e 2 ^ 3 r 5 ^ 3 /vq, ex- 
cept for r 3> Vat. Since r ~ vot is nearly the maximum 
particle separation that can be achieved in the time t, 
only a very tiny large-r tail will experience a different 
eddy-diffusivity than this. In Fig. [7] we also compare 
the Monte Carlo results for the KS diffusion model with 
the exact parameter- free predictions Q117[) . l|118[l of the 
power-law diffusion model (|114[) for d = 3 and ( = 5/3. 
The agreement is reasonably good. Furthermore, the 
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FIG. 6: Numerical results for (\r(t) - r | 2 > in the KS diffu- 
sion model (|128|l . fczv/fci = 10 21 : Monte Carlo solution of the 
diffusion equation (o) with Ca = 1, S = 10 4 and mean-field 
(— ). Also MC results (• ■ ■ ) for inertial-range model ([84]l . 




FIG. 7: Logarithm of the rescaled PDF of pair-separations 
at 23 different times in the t 6 range, for KS diffusion model 
with kn/ki = 10 21 . Monte Carlo results for S = 10 6 and 
Ca = 1. The straight line marked with circles (o) is the 
analytical result (|116[) for d — 3 and 7 = 1/3. 



Monte Carlo results for the KS diffusion model agree 
almost perfectly with the Monte Carlo results for the 
power-diffusion model presented in section IIV Al This 
suggests that if our Monte Carlo could be carried out 
with a small, constant time-step in a long t -range, then 
the PDF would approach the exact self-similar form of 
the power-law diffusion model. 

We have conjectured that the growth laws of the KS 
models themselves, asymptotically for fcjv/fci 3> 1, are 
the t 2 , t 1 , t 6 and t 1 powers that we have found in the KS 
diffusion models. How can this be reconciled with the t 9 / 2 
law predicted in [H and verified to greater or lesser extent 
in subsequent KS simulations [1, [13, [ll|? We argue that 
the observed t 9 ^ 2 is an artefact of the modest fcjv / k\ ra- 
tios achieved in these simulations, which tends to "blend" 
the distinct scaling ranges, in particular the early-time 
t 1 and t 6 ranges, between which lies a broad transition 
zone. In support of this argument, we have performed 
a sequence of Monte Carlo simulations of the KS diffu- 
sion model with scale ratios fcjv/fci = 10 3 , 10 4 , 10 5 , 10 6 . 
The last ratio is chosen to correspond roughly to that em- 
ployed in the previous KS simulations |slnj,lll|]. Because 
the range of time-scales is not so great, we have been able 
to carry out the time-integration not only with the adap- 
tive algorithm employed up until now, but also with a 
constant time step At = 0.1— which resolves the effects 
of even the smallest eddies, equivalent to that used in re- 
cent KS simulations [HI EH ■ The results of the two time- 
advancement schemes for the dispersion curves are identi- 
cal when plotted in log- log. As illustrated in Fig. [SI a t 9 / 2 
regime seems to appear as we increase the ratio fcw/fci. 
This figure should be compared with Fig. 2 of [l(| and 
Fig. 1 of 11], which it matches very closely. Although 
we see a similar "t 9 / 2 -range" at the values of kjq / k\ used 
in previous KS simulations, covering 1-2 decades in time, 
it is clear from our results in Fig. [S] that this is only a 
transitional regime of the KS diffusion model. In fact, for 
the case fcjv/fci = 10 6 which shows the long "i 9 ' 2 -range" 
we find t to = 10 -2 and thus the broad transition zone 
between the t 1 and t 6 laws covers the interval from 10~ 4 
to 10°. This includes all of the apparent "i 9 / 2 -range" . If 
we go to fcjv/fci = 10 8 , the power-law steepens into a 
t 5 -law. At still larger values of fcjv/fci four asymptotic 
scaling ranges emerge, with distinct power-law scalings 
of t 2 , t 1 , t 6 and t 1 . We expect that the same is true of the 
KS models themselves at sufficiently large kx/k\. 

Finally, we note that for k^/ki < 10 4 , the short 
range of supcrdiffusive growth of dispersion approxi- 
mates a £ 3 -law. This agrees with the observations of 
@, HH f° r models. Note, however, that the physics is 
completely different from turbulent Richardson diffusion, 
which would allow t 3 ranges of arbitrary extent. In fact, 
the narrow range of such a power-law in our KS diffu- 
sion model arises only because of the "merging" of many 
distinct ranges. In particular, the exponent of the ap- 
parent power-law must decrease with decreasing k^/ki 
to match the ^-law starting at r = L, until finally the 
supcrdiffusive range disappears entirely when kj^/ki sa 1. 
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Q for ensembles with large mean velocities. 

Synthetic models of turbulence such as Kinematic Sim- 
ulations have been used to investigate turbulent trans- 
port of passive objects (particles, lines, etc.) in such var- 
ied problems as environmental flow, aeroacoustics, kine- 
matic magnetic dynamo, and supcrfluids [? ? ]. How- 
ever, such numerical studies must clearly be employed 
with utmost caution, especially to derive conclusions 
about turbulent transport at very high Reynolds num- 
bers. The difference in sweeping effects in synthetic Eu- 
lerian turbulence and in real hydrodynamic turbulence 
imply not only quantitatively different scaling laws but 
also substantially different physics. 
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Appendix A: Monte Carlo Time-Step 



FIG. 8: Monte Carlo results for (|r(t)| 2 ) in KS diffusion model 
for various values of fav/fci, with At = O.l^j-, S = 10 4 . 



We tested the dependence of the log-log plots of dis- 
persion on the value of Ca. Wc plot in Fig. |2]the Monte 
Carlo results for values of Ca ranging from 1 to 1CP 6 . 
There is no observable change in the behavior. 



V. CONCLUSIONS 

We have derived in this paper a diffusion equation for 
particle-pair dispersion in synthetic Eulerian turbulence 
modelled by Gaussian velocity ensembles. The main an- 
alytical result is the formula (|54j) for the 2-particle dif- 
fusivity and its special cases (|66[) for frozen velocities 
and (|76|) for finite time-correlated velocities. Although 
the description of pair-dispersion as a diffusion process 
is not exact (except in certain limiting cases), it arises 
from a well-motivated set of analytical approximations. 
Our results confirm the physical argument of Thomson 
& Devenish Q that pair-dispersion in such models is 
fundamentally altered by sweeping dccorrelation effects, 
not experienced by particle pairs in hydrodynamic turbu- 
lence. Thus, the t 3 -law observed in previous simulations 
with synthetic turbulence [H-Q is quite likely an arte- 
fact either of the numerical approximations employed or 
of the shortness of the inertial ranges. However, we ar- 
gue as well for a similar origin of the i 9 / 2 -law proposed 
by Thomson & Devenish Q for synthetic turbulence en- 
sembles with zero mean velocities. Solutions of our diffu- 
sion model for such ensembles at Reynolds numbers com- 
parable to those employed in KS simulations that show 
a £ 9 / 2 -law range reproduce that finding, but our model 
yields instead distinct t 2 , t 1 , t 6 and ^-ranges at higher 
Reynolds numbers. We thus argue that the asymptotic 
high Reynolds-number behavior of particle dispersion in 
synthetic Eulerian turbulence with zero mean-velocities 
is the same as that predicted by Thomson & Devenish 



Appendix B: Number of Fourier Modes 

We also tested the dependence of our dispersion results 
for the KS diffusion models on the number of Fourier 
modes N. We show in Fig. [10] log-log plots of the dis- 
persion curves for different values of N, obtained from 
Monte Carlo calculations with Ca = 1 and S = 10 4 . 
The results are nearly indistinguishable for N > 100. All 
of our simulations in the text used TV = 500. 
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FIG. 9: Monte Carlo results for (|r(i) - r | 2 ) in the inertial- 
range diffusion model calculated with S = 10 2 samples and 
varying Ca = 1 to Ca = 10~ 6 . 



FIG. 10: Monte Carlo results for (\r{t) - r | 2 ) in the KS 
diffusion model calculated with Ca = 1, S = 10 4 samples, 
varying number of Fourier modes from N = 10 to N = 10 4 . 
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[31] This is the first point where we have invoked incompress- 
ibility. 

[32] Here the superscript in D " is used to indicate the spa- 
tial scaling exponent a for which the constant in (J3]) is 
calculated. Using eqs. (2.14) and (2.16) in p2] to calcu- 
late D\ gives 

£ q r(^)r(^f±2) 

°° V 8 a + 1 T (2=^) r (^±f±2) 



[33] The equations of [? ] are unfortunately marred by several 
misprints 

[34] The standard formula Cl = j^r, — , 2 \ ■ ? — t^Ck re- 

lates the constants Cl, Cjr for a fc~' 1+a ' power-law spec- 
trum; e.g. see Q, eq.(13.100). This leads to Cl = 1.9727 
and C' L — ClJoo = 1.262, the choice of the previous two 
subsections. 



